/*  sdfcbrt.cpp  by K.Tsuru  */
// function ID 3010     DRADIX only
/***************************
SDouble class
cubic root of a
It returns a*RecCbrt(a*a).
****************************/
#ifndef SN_H
#include "sn.h"
#endif
static const char* const func = "Cbrt";

SDouble Cbrt(const SDouble& a){
    SDouble d, y;

//When 'a' is a short number including |a| = 1.0 or 0.0
//It costs a little time bacause the "db" is a short number.
    uint xfig = a.Last() - a.First();
    if(xfig < 4){
        SDouble db;
        y.SNError();    //reset error
        db = pow( fabs(doubleD(a, 0)), 1.0/3.0);
        if(y.SNError()==y.NO_ERR){  //detect an error in doubleD()
            if(a.Sign(3010)<0) db.ChangeSign();
            y = a - db*(db*db);
            if(y.Sign(3007)== 0) return db;
        }
    }
    y = a*a;
    y = a*RecCbrt(y);
/**************************************************
It adds a correction. Maybe this is not necessary.
Let d = a - y'^3
 y = a^(1/3) =(y'^3+d)^(1/3) = y'{1+d/(y'^3)}^(1/3)
 = y' + d/{3*(y'^2)}+...
****************************************************/
#if 0
    if(y.PreferSpeed() == OFF){
        SDouble y2;
        y2 = y*y;   d = a - y2*y;
        if(d.Sign()){
            d = d/y2;
            y += DsDiv(d, 3);
        }
    }
#endif
    if(y.Verify()){     // check |a - y^3| << a ?
        d = a - (y*y)*y;
        if(!d.IsMLT(a)){  //It sees the relative error. d << a ?
            y.SetError(y.VERIFY, func, 3010);
        }
    }
    return y;
}